Settlement on asphalt concrete pavement in permafrost regions under the dynamic load of aircraft

The settlement in permafrost regions has a significant effect on the safety of the aircraft. Therefore, a numerical model of asphalt concrete pavement and subgrade in permafrost regions is proposed and verified by comparing with previous studies. Numerical models under the dynamic load of aircraft are used to analyze the settlement nephogram, time-dependent curve, and settlement curve. Results show that the influences of different parameters on settlements decrease slowly at the depth of 0–1.45m, then decrease quickly at the depth of 1.45-2m, finally tend to be stable. The peaks of settlements increase with the growth of international roughness index, taxiing speed, and aircraft’s weight. The settlement increases little when the international roughness index is more than 3. The settlement varies significantly when the taxiing speed is from 30m/s to 70m/s. The study provides a theoretical basis for the construction and maintenance of asphalt concrete pavement and subgrade in permafrost.


Introduction
A modern international corridor of the airport system oriented to "Belt and Road" will be gradually built to promote the connectivity between China and its surrounding countries in the future [1]. Hence, 68 new and renovated general airports will be built in permafrost regions in accordance with the medium-term and long-term plan for airport construction in China [1]. However, the permafrost is thermally unstable, and is easy to cause the thawing and settlement with the effect of pavement weights and aircraft loads [2]. Moreover, the investigation results show that engineering structures in permafrost regions suffer from cracks, deformation, collapse caused by settlement. Therefore, it is urgent to settle the matter of settlement in the construction of airport in a permafrost region affected by the temperature, pavement weights, and aircraft loads.
Scholars had studied the settlement of permafrost [3,4]. Morgenstern and Nixon [5], and Foriero and Ladanyi [6] applied Terzaghi consolidation theory, the large strain thawing consolidation theory, and heat conduction equation to analyze the settlement of the saturated and unsaturated soil. Hou [7] adopted conservation law of momentum and mass, entropy inequality, and the mixture theory to established the settlement model. They considered the influence of the volume fraction condition. With the development of computer technology, Fortier et al. [8], and Wang [9] applied the water-thermal-mechanical coupling model to analyze the influence of subgrade parameters, snow cover, slope orientation on the settlement. Their study showed that the thawing inter-layer, slope effect, and thawing depth had the significant impact on the uneven settlement. Liu [10], and Zhao [11] explored the influence of subgrade parameters and freeze-thaw cycles on the mechanism of railway settlement based on the water-heat-mechanical response model. Vladimir and Elena [12] had successfully predicted the response of settlement in the Russian northern ports based on the numerical simulation method. In addition, scholars had also carried out a large number of experimental studies. Zaretskii et al. [13], Arenson et al. [14], Zheng [15], and Kenan [16] analyzed the connection of settlement, dry density, moisture content, and freezing temperature based on the triaxial compression test, volume change test, unconfined compression test, and freeze-thaw cycle test. Peng et al. [17], Zhu [18], and Lu et al. [19] used the in-situ data of Qinghai-Tibet highway to conclude four types of settlement patterns, including low steady subsidence, fluctuating subsidence, step subsidence and high steady subsidence. Their research explained that the wave of pavement, and the longitudinal crack and collapse of subgrade mainly occurred in the first freeze-thaw cycle stage and were caused by the longitudinal and lateral settlement. Ma [20] and Chen [21] studied the evolution law of settlement using the in-situ monitoring data from the Qinghai-Tibet Highway test. Wang [22] explored the mechanism of thawing settlement of tunnel in a permafrost region, and discovered that the thawing settlement primarily took place in the early stage. Dong [23] and Wang et al. [24] analyzed the mechanism of settlement of pipeline by the field trial. The above research focuses on railway, highway, tunnel, and pipeline engineering, and does not fully considered the characteristics of pavement in the airport including wide width, short length, low pavement roughness, and large aircraft's weight. Meanwhile, the existing research on the mechanism of settlement has more consideration on static factors including subgrade parameters, slope, and freeze-thaw cycle, but less consideration on the dynamic load of aircraft. However, the settlement of asphalt concrete pavement in permafrost regions is affected by multi-dimensional factors, such as self-weight stress from subgrade and pavement structures, subsidiary stress from dynamic load of aircraft, and temperature field from environmental factors. Moreover, the dynamic load of aircraft is affected by the aircraft's weight, international roughness index, and taxiing speed. Therefore, based on the characteristics of airport and the actual conditions of permafrost subgrade in China, this paper establishes aircraft and pavement response model to reveal the settlement mechanism of permafrost subgrade under the dynamic load of aircraft. The study provides a theoretical basis for the construction and maintenance of asphalt concrete pavement in a permafrost region.

Heat conduction equations
Because the asphalt concrete pavement is disturbed by solar precipitation and radiation, water content of structural layer and geological lithology, some assumptions are made to ensure computational efficiency and effectiveness. Firstly, the pavement, natural surface, and subgrade are uniform and isotropic. Secondly, the boundary of pavement, natural surface, and subgrade does not consider the migration of moisture. Thirdly, the layers of pavement, natural surface, and subgrade meet the law of energy conservation, and only include the effects of heat transfer. Therefore, the three-dimensional unsteady heat conduction express as Eq (1) during freezing and thawing [25,26]: where C and T are the heat capacity and the temperature, respectively. t and L are the time and the transformed latent heat of water, respectively. ρ i and w i are the density and volume fraction of ice, respectively. k x , k y , and k z are the components of thermal conductivity. Meanwhile, the three-dimensional mass transfer express as Eq (2) with evaporation and other factors being ignored [26,27]: where w u and ρ w are the volume fraction of unfrozen water and the density of water, respectively. K x , K y , and K z are the components of hydraulic conductivity. φ is the total potential energy, and equals to ϕ+h. ϕ and h are the volumetric potential energy and the gravitational potential, respectively. h can be ignored. The content of unfrozen water can be shown in Eq (3) [26]: Eqs (1)-(3) can applied to obtain Eq (4) [26]: The volumetric potential of unfrozen water is shown in Eq (5) [26]: Eq (6) can be obtained with the differential water capacity and water diffusion coefficient into Eq (5) [26]. Based on Eqs (4) and (6), the three-dimensional unsteady heat conduction express as Eq (7) [26].
where D and K are the coefficient of water diffusion and the hydraulic conductivity, respectively.

Boundary conditions
Based on the surface temperature in a permafrost region, the first boundary condition is established in Eq (8) and is applied on the surface of numerical model [21,26]: where T 0 is the mean annual ground temperature of the boundary layer, taken as the mean annual air temperature, ΔT is the increment in the temperature of the boundary layer, and is set to 2.5˚C for the natural surface and 4.5˚C for the pavement. a is the amplitude of variation in temperature, and is set to 12˚C for the natural surfaces and 15˚C for the pavement. t is the time in hours, b is the initial phase, equals to zero at the initial set time of April 1 in northeast China. α is the increment in the global temperature over the next 30 years, and is set to 1.5˚C. The second boundary condition is an adiabatic boundary condition in Eq (9), and is applied on both sides of numerical models [16,26]: The third boundary condition is a constant temperature gradient of 0.03˚C/m. This is applied at the bottom of the numerical model, as shown in Eq (10) [26,28]: The support of landing gear is completely rigid, and the axial motion has a linear damping. (7) The dual wheels on one side are equivalent to one equivalent single wheel. Based on the theory of time-varying mechanics, the aircraft and pavement are regarded as two independent subsystems. The contact parts between wheel and pavement are connected by the deflection coordination condition and the force balance condition. According to the vibration theory, the vibration equation of aircraft is established by combining the Lagrange equation with the force balance, moment balance and displacement coordination of system, as shown in Eqs (11)-(25) [29].

Aircraft vibration equations
where K E is the kinetic energy of the aircraft system. P E is the potential energy of the aircraft system. D E is the energy consumption of the aircraft system. M 0 is the mass of the aircraft. M 1 , M 2 , and M 3 are the mass of left and right main landing gear and front landing gear, respectively. Z 0 is the vertical displacement of aircraft affected by the international roughness index. Z 1 , Z 2 , and Z 3 are the vertical displacement of left and right main landing gear and front landing gear, respectively. J x and J y are the mass moment of inertia of X and Y axis, respectively. θ x and θ y are the rolling rotation and pitching rotation, respectively. K xHZ and K xQZ are the suspension stiffness of left and right main landing gear and front landing gear, respectively. C xHZ and C xQZ are the suspension damping of left and right main landing gear and front landing gear, respectively. K LHZ and K LQZ are the tire stiffness of left and right main landing gear and front landing gear, respectively. C LHZ and C LQZ are the tire damping of left and right main landing gear and front landing gear, respectively. Z 4 , Z 5 , and Z 6 are the tire vertical displacement of left and right main landing gear and front landing gear, respectively. Z 7 , Z 8 , and Z 9 are the pavement roughness corresponding to left and right main landing gear and front landing gear, respectively. L H and L Q are the distance between center of gravity and the rear and front axle, respectively. L y is the distance between left and right axle.

Establishing the numerical model
In Fig 2, the natural surface was made of 0.4-m-thickness clay soil, 1.6-m-thickness silty clay soil, and 18-m-thickness strongly weathered rock based on an engineering geological survey [1,26]. The specifications of the design of asphalt pavements for civil airports indicated that the pavement structure was made of upper and under surface, upper and under base layer, and subbase layer [26]. To eliminate the boundary effect, a 40-m-wide natural surface was placed on both sides of asphalt concrete pavements [26]. The numerical model of asphalt concrete pavement was 45-m-width and 15-m-length [26]. The depth was along the Y direction, width was along the X direction, and length was along the Z direction [26]. Eqs (8), (9) and (10) were applied to the boundary conditions of numerical models. The bottom of strongly weathered rock was constrained along the depth, length and width. The edges of strongly weathered rock, silty clay soil, and clay soil were denied the freedom degree of length and width direction [26]. The friction coefficient of interface layer was 0.5. The calculation time was 10 years [26]. The initial temperature field was -1˚C [26,30]. The phase change in a certain range of temperature was simulated by the apparent heat capacity method [26]. Tables 1 and 2 showed the thermal parameters of weathered rock, silty clay soil, clay soil, and pavement [1,26]. The elastoplastic behavior of the soils was described by the Mohr-Coulomb constitutive model. Based on the research of Li [1] and Liu et al. [26], the mechanical parameters of permafrost were related to the temperature, and solved by Eqs (26)- (29). Table 3 showed the coefficients of expansion of strongly weathered rock, silty clay soil, and clay soil [1,26]. Tables 4 and 5 showed the mechanical parameters of the pavement structure [1,26].
The expansion coefficient is a dimensionless index. The three-dimensional mass-spring-damper system with six degrees of freedom (DOFs) was applied to simulate the aircraft, as shown in Fig 1 and Table 6 [31]. The body of aircraft was three DOFs, such as the vertical displacement at the mass center, rolling rotation and pitching rotation. Each tire was a lumped mass with vertical displacement. The aircraft was described by beam elements (MPC184). The suspension and nonsuspension systems and the inertia moments of pitching and rolling rotations were described by mass elements (Mass21). The spring and damper were described by a spring-damper element (Combine14). The dynamic load of B737-800 aircraft was described as a uniformly distributed rectangular load, and arranged symmetrically along the centerline of the upper surface. The elements of the loading area were selected at 0.30-m-width and 0.43-m-length according to the widths of tire ribs and grooves [32]. The tire pressure was calculated by the loading area and tire load. The moving load was described by controlling the loading area on pavement surface and the specific interval based on the speed. The speed of aircraft was from 10m/s to 70m/s in the numerical model.

Verifying the numerical model
Owing to a lack of testing data on runways, the reliability of numerical model was verified by a comparison with the temperature field, frozen layer depth, and settlement in previously studies. A transient thermal analysis of runway was conducted over 10 years with the initial temperature field. Figs 3 and 4 showed the temperature field and the depth of frozen layer after 10 years of operation [26]. Fig 5 showed the displacement field after one year of operation [26]. Fig 3 showed that the maximum depth of the frozen layer was on the 351 th day, at the end of March in the second year, which were in agreement with the study of Liu et al. [33], Wang [34] and Liu et al. [26]. Fig 4 showed that the pavement started to freeze in the middle of October, reached the maximum depth at the end of March, began to thaw at the beginning of April, and had completely thawed by the end of April, which were consistent with the report from Li [1], Duan [30] and Liu et al. [26]. Moreover, the maximum depth of the frozen layer was about 3.35m in runway engineering and deeper than that of highway and railway engineering, which The maximum settlements at the center of the pavement was -0.023m, and the maximum settlements at the center of the natural surface was -0.091m in runway engineering, which were in accord with the result from LeBlanc et al. [37] and Liu et al. [26]. Hence, the parameters and boundary conditions of runway models were reasonably applied to analyze the settlement response of runway in a permafrost region.

Evolution of settlement
To observe the evolution of settlement, the temperature field nephograms on the first day of each month in the first year were extracted in Fig 6. In Fig 6(A), January was the coldest month and the surface of the pavement was significantly influenced by the external environment and began to cool down firstly. Hence, the temperature of the pavement increased with the depth increased, and was barely effected by the external environment when the depths were more than 12m. In Fig 6(B), the air temperature increased from February to April. In April, the temperature of the pavement center exhibited high-low-high-low characteristics with the increase of the depth. The surface of the pavement began to warm up firstly with the growth of the air temperature. The temperature of the  pavement changed more and more slowly with the increase of the depth, and the high-temperature interlayer became smaller gradually in the pavement structure. In Fig 6(C), the air temperature continued to increase from May to July. As July was the hottest month, the temperature of the pavement center exhibited high-low characteristics with the growth of the depth. The high-temperature interlayer disappeared as the whole pavement was in high-temperature. In Fig 6(D), the air temperature gradually decreased from August to October. In October, the temperature of the pavement center started to slow and exhibited low-high-low characteristics with the increase of the depth. The surface of the pavement was effected by the external environment firstly, and began to slow in temperature. The temperature of the pavement became slower as the depth increased, and the high-temperature interlayer could be clearly seen in the pavement structure. In Fig 6(E), because the air temperature continued to decrease from November to December, the cooling depth of the pavement increased and the high-temperature interlayer was deepening and expanding.

Settlement under different loads
To observe the influence of different loads on the settlement, the settlements under different loads along X-axis and depth at the middle of August of the tenth year were shown in Fig 7. In Fig 7, the dynamic load of aircraft was obtained at the international roughness index of 3, the taxiing speed of 30m/s, and the aircraft's weight of 57878kg. In Fig 7(A), the settlement of self-weight of pavement was the smallest. The settlement of dynamic load was larger than that of static load. Because the settlement of permafrost was mainly influenced by temperature and load, which was a coupling process of temperature and load. When the warm season came, the frozen ice in the soil thawed into water to reduce the volume and strength of soil. The permafrost was subject to settlement because of the deadweight of the pavement structure and subgrade. Meanwhile, the load of aircraft also caused the settlement of the permafrost. Moreover, because the dynamic load included the self-weight of aircraft and the additional load caused by the international roughness index and taxiing speed, the dynamic load of aircraft was greater than that of static load. And the settlement of dynamic load was larger than that of static load. In Fig 7(B), the shape of settlement along the depth was similar at different loads. Firstly, the settlement decreased slowly at the depth of 0m to 1.45m, then decreased quickly at the depth of 1.45m to 2m, finally tended to be stable at the depth of 2m to 11m. When the depth was from 0m to 1.45m, the settlement of dynamic load was the largest and the settlement of static load was bigger than that of self-weight of pavement. When the depth was from 2m to 11m, the difference of settlement was little caused by dynamic load, static load, and self-weight of pavement. Because the elastic modulus of permafrost subgrade varied greatly because the temperature caused the settlement in a relatively shallow range. Meanwhile, the stress induced by the load decreased gradually along the depth direction, so the settlement caused by the pavement structure and aircraft also reduced gradually along the depth direction.

Influence of different parameters on the settlement
In the aircraft and pavement response analysis of numerical models, the international roughness index of 1 was to represent a new pavement or the best pavement condition, the international roughness index of 3 was to represent the good pavement condition with a little damage, and the international roughness index of 5 was to describe the poor pavement condition with much damage [29]. The speed of 10m/s indicated that the aircraft was beginning to taxiing, the speed of 30m/s indicated an accelerated taxiing of aircraft, and the speed of 70m/s indicated that the plane was about to take off. The weight of 42275kg was the weight of aircraft at no load, the weight of 57878kg was the weight of aircraft at some load, and the weight of 73482kg was the weight of aircraft at full load. Fig 8 represented the time history curve of settlement at different positions of pavement when the taxiing speed was 30m/s and the weight was 57878kg. In Fig 8(A), when the international roughness index was from 1 to 5, the shapes of time history curves were similar with different peaks of time history curves. When the time was from 1 st year to 5 th year, the settlement had a rapid increase. While the time was from 5 st year to 10 th year, the settlement had a minimal change. The reason might be that the consolidation settlement of soil mainly occurred in the previous years, and the subsequent consolidation settlement gradually decreased. Meanwhile, the voids and moisture in the soil decreased gradually to reduce the settlement of soil. When the international roughness index varied from 1 to 3 at the same time, the settlement increased quickly. While the international roughness index varied from 3 to 5 at the same time, the settlement grew little. It was because that the dynamic load of aircraft caused by damaged pavement was much greater than that caused by good pavement. In Fig 8(B), the development laws of time history curves were similar to those of Fig 8(A), and the peaks of time history curves were smaller than those of Fig 8(A). The explanation of the phenomenon was that the load at the center of the wheel track was greater than that at the center of the pavement, and the settlement at the center of the wheel track was larger than that at the center of the pavement.
The top of the upper surface layer was set to 0m. Fig 9(A) showed the influence of international roughness index and years on the settlement at the depth of 0m-11m. Fig 9(B) showed the influence of international roughness index and years on the settlement at the depth of 2m-11m. Fig 9(B) was an enlarged view of Fig 9(A) at the depth of 2m-11m. As the time and the international roughness index varied, the change law of settlement was similar along the depth. Firstly, the settlement decreased little with the depth being from 0m to 1.45m, then changed large with the depth being from 1.45m to 2m, and finally tended to fluctuate between 0.0055m and 0.0125m with the depth being from 2m to 11m. The phenomenon showed that the settlement mainly appeared in the soil instead of the pavement structure layers, and the significant area of settlement was below the pavement structure layers. Therefore, protecting the soil temperature stability below the pavement structure layers was helpful to control the settlement of runway. Fig 10 described the time history curve of settlement at different positions of pavement when the international roughness index was 3 and the weight was 57878kg. In Fig 10(A), the change laws of time history curves were similar with different peaks of time history curves, which were different with the increase of the taxiing speed. Because the reduction of voids and moisture in the soil mainly occurred in the previous years and the subsequent reduction was small and slow, the settlement after 5 years of operation in the northeast of China varied little at different taxiing speed. Moreover, because the taxiing speed was beneficial to increase the dynamic load of aircraft, the settlement grew with the increase of the taxiing speed. In Fig 10  (B), the shape of time history curves were the same to those of Fig 10(A), and the peaks of time history curves were smaller than those of Fig 10(A). The reason might be that the settlement caused by the wheel at the center of the wheel track was greater than that at the center of the pavement. The varied value of settlement at the taxiing speed of 10m/s to 30m/s was significantly smaller than that at the taxiing speed of 30m/s to 70m/s. The phenomenon indicated that the faster the taxiing speed, the larger the settlement caused by the dynamic load of aircraft. Fig 11 showed the influence of taxiing speed on the settlement at different years and depth. As the time and the taxiing speed varied, the change law of settlement was similar along the depth. The settlement had little change with the depth being from 0m to 1.45m, varied largely with the depth being from 1.45m to 2m, and kept stable in a small range with the depth being from 2m to 11m. Because the thickness of the pavement structure layer was 1.45m, and the settlement mainly appeared in the subgrade and decreased as the depth increased. Hence, the temperature stability of soil in a certain range below the pavement structure was the key to control the settlement of runway. Fig 12 represented the time history curve of settlement at different positions of pavement when international roughness index was 3 and the taxiing speed was 30m/s. In Fig 12(A), time history curves of settlements had the similar shapes and different peaks with the growth of aircraft's weight. The settlement was positively correlated with the aircraft's weight. When the aircraft's weight varied from 42275kg to 73482kg, the settlement tended to be stable after 5 years of operation in the northeast of China. Moreover, the settlement was more than 3cm at the aircraft's weight of 73482kg after 5 years of operation in the northeast of China. In Fig 12(B), the laws of time history curves were similar and the peaks of time history curves were different with the aircraft's weight being from 42275kg to 73482kg, which was similar to those of Fig 12  (A). However, the settlement was smaller than that of Fig 12(A) at the same time and aircraft's weight. The reason was that the load of aircraft decreased with the growth of distance away from the wheel track, and the settlement was positively correlated with the aircraft's weight. Fig 13 showed the influence of aircraft's weight on the settlement at different years and depth. when the time and the aircraft's weight increased, the settlement decreased slowly with the depth being from 0m to 1.45m, reduced quickly with the depth being from 1.45m to 2m, and tended to be the range of 0.0055m to 0.0125m with the depth being from 2m to 11m. Because the material of pavement structure layer had less voids and moisture than that of subgrade, and the void and moisture were helpful to cause the settlement of runway. Hence, controlling the void and moisture was also beneficial to keep the stability of subgrade in a permafrost region.

Conclusions
This study proposed and verified a numerical model of asphalt concrete pavement and subgrade in permafrost regions by the temperature field, frozen layer depth, and settlement comparison of previous research without considering the dynamic load of aircraft. This study also analyzed the settlement under different parameters. The main conclusions were shown as follows: 1. In January, the temperature of pavement decreased as the depth increased, and was barely effected by the external environment when the depths were more than 12m. In April, the temperature of pavement center exhibited high-low characteristics with the growth of the depth. The surface of pavement began to warm firstly as the air temperature increased. The temperature of pavement slowed with the increasing depth, and an obvious low-temperature interlayer appeared in the pavement structure. In July, the temperature of pavement center exhibited high-low characteristics with the depth increased. The low-temperature interlayer disappeared. In October, the temperature of pavement center started to slow and exhibited low-high-low characteristics with the depth increased. The temperature of pavement became slower with the increase of the depth, and the high-temperature interlayer could be seen in the pavement. When it was from November to December, the cooling depth of pavement increased and the high-temperature interlayer was deepening and expanding.
2. The shapes of settlements under self-weight of pavement, static load, and dynamic load were similar and reduced gradually along the depth direction. The settlement of static load was larger than that of self-weight of pavement, and smaller than that of dynamic load. The settlements of different loads decreased slowly at the depth of 0m to 1.45m, then decreased quickly at the depth of 1.45m to 2m, finally tended to be stable at the depth of 2m to 11m.
3. When the time was from 1st year to 5th year, the settlement had a rapid increase. While the time was from fifth year to tenth year, the settlement had a minimal change. The settlement at the center of the wheel track was larger than that at the pavement center. The laws of time history curves were similar with different peaks of time history curves for different international roughness index, taxiing speed, and aircraft's weight. As the time varied, the change law of settlement was similar along the depth for different international roughness index, taxiing speed, and aircraft's weight. Firstly, the settlement decreased little with the depth being from 0m to 1.45m, then changed large with the depth being from 1.45m to 2m, and finally tended to fluctuate between 0.0055m and 0.0125m with the depth being from 2m to 11m.
4. When the international roughness index varied from 1 to 3 at the same time, the settlement increased quickly. While the international roughness index varied from 3 to 5 at the same time, the settlement grew little. The varied value of settlement at the taxiing speed of 10m/s to 30m/s was significantly smaller than that at the taxiing speed of 30m/s to 70m/s. When the aircraft's weight varied from 42275kg to 73482kg, the settlement increased in the first five years and tended to be stable after 5 years of operation in the northeast of China. Moreover, the settlement was more than 3cm at the aircraft's weight of 73482kg after 5 years of operation in the northeast of China.